dir_data <-"../Data/lidar"
year <- "2013"
dirname <- paste(dir_data,"vector",year,"treetops",sep="/")
ws <- 18
test_filename <- paste0("treetops_lmf_ws",ws)
test_filename_full <- paste(dirname,test_filename,sep = "/")
# filename <- "../Data/lidar/vector/2013/treetops/treetops_lmf_ws18"
dir_crowns <- paste(dir_data,"raster","crowns",sep="/")
dir.create(dir_crowns,recursive = T)
## Warning in dir.create(dir_crowns, recursive = T): '../Data/lidar/raster/crowns'
## already exists
## Read underlying raster
raster_filename <- paste(dir_data,"raster","raster2013.tif",sep="/")
r <- raster(raster_filename)
treetops <- readOGR(test_filename_full)
## OGR data source with driver: GeoJSON
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/2013/treetops/treetops_lmf_ws18", layer: "2013_ws18"
## with 2510 features
## It has 2 fields
plot(r)
plot(treetops,add=T)

filename_dalponte <- paste(dir_crowns,
paste0(
paste(year,ws,"dalponte",
sep = "_"),
".tif"
),
sep="/")
# crowns_dalponte <- dalponte2016(
# r,
# treetops,
# th_tree = 2,
# th_seed = 0.45,
# th_cr = 0.55,
# max_cr = 50,
# ID = "treeID"
# )()
# writeRaster(crowns_dalponte,filename_dalponte)
crowns_dalponte <- raster(filename_dalponte)
plot(crowns_dalponte,col=pastel.colors(1000))

filename_silva <- paste(dir_crowns,
paste0(
paste(year,ws,"silva",
sep = "_"),
".tif"
),
sep="/")
# crowns_silva <- silva2016(r,treetops)()
# writeRaster(crowns_silva,filename_silva)
crowns_silva <- raster(filename_silva)
plot(crowns_silva,col=pastel.colors(1000))

# contours <- rasterToPolygons(crowns_dalponte,dissolve = TRUE)
subsampling <- 100000
pixel_area <- 0.5*0.5
ncells <- r@ncols*r@nrows
hist <- hist(crowns_dalponte,breaks=2510,plot=F, maxpixels=subsampling)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.
areas <- sort(hist$counts)*ncells/subsampling*pixel_area
plot(areas,pch=".")
hist <- hist(crowns_silva,breaks=2510,plot=F, maxpixels=subsampling)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.
areas <- sort(hist$counts)*ncells/subsampling*pixel_area
points(areas,pch=".",add=T,col="red")
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter
